R programming: generics

M1 MIDS/MFA/LOGOS

Université Paris Cité

Année 2024

Course Homepage

Moodle

stopifnot(
  require(Hmisc),
  require(skimr),
  require(patchwork),
  require(ggforce),
  require(glue),
  require(ggfortify),
  require(broom),
  require(tidyverse)
)

tidymodels::tidymodels_prefer(quiet = TRUE)

old_theme <-theme_set(
  theme_minimal(base_size=9, 
                base_family = "Helvetica")
)

gc <- options(ggplot2.discrete.colour="viridis")
gc <- options(ggplot2.discrete.fill="viridis")
gc <- options(ggplot2.continuous.fill="viridis")
gc <- options(ggplot2.continuous.colour="viridis")

Objectives

stopifnot(
  require(rlang),
  require(lobstr),
  require(sloop),
  require(devtools),
  require(usethis),
  require(testthat),
  require(generics)
)

Loading required package: rlang

Loading required package: lobstr

Loading required package: sloop

Loading required package: devtools

Loading required package: usethis

Loading required package: testthat

Loading required package: generics

Generics and S3 classes

Let us first create an instance of class lm.

lm0  <- lm(Gas ~ Insul * Temp, MASS::whiteside)

Question

• What does function class() do?
• Is it possible to belong to type list and to class lm simultaneously?
• In R what is an attribute?
• How do we set and get attributes?
  
• What does function inherits() do?

Solution

class(lm0)

[1] "lm"

is.list(lm0)

[1] TRUE

attributes(lm0)

$names
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "contrasts"     "xlevels"       "call"          "terms"        
[13] "model"        

$class
[1] "lm"

inherits(lm0,"list")

[1] FALSE

inherits(lm0, "lm")

[1] TRUE

Question

Load package sloop.

• What does sloop::otype() do? Apply it to an object of class lm.
• What happens when we first unclass() the object?

Solution

sloop::otype(lm0)

[1] "S3"

sloop::otype(unclass(lm0))

[1] "base"

Question

sloop exports functions s3_class() and s3_get_method()

• Apply s3_class() to all members of lm0
• What is the otype of autoplot() applied to an object of class lm?

Solution

sloop::s3_class(lm0)

[1] "lm"

sloop::s3_get_method(autoplot.lm)

function (object, which = c(1:3, 5), data = NULL, colour = "#444444", 
    size = NULL, linetype = NULL, alpha = NULL, fill = NULL, 
    shape = NULL, label = TRUE, label.label = ".label", label.colour = "#000000", 
    label.alpha = NULL, label.size = NULL, label.angle = NULL, 
    label.family = NULL, label.fontface = NULL, label.lineheight = NULL, 
    label.hjust = NULL, label.vjust = NULL, label.repel = FALSE, 
    label.n = 3, smooth.colour = "#0000FF", smooth.linetype = "solid", 
    ad.colour = "#888888", ad.linetype = "dashed", ad.size = 0.2, 
    nrow = NULL, ncol = NULL, ...) 
{
    p1 <- p2 <- p3 <- p4 <- p5 <- p6 <- NULL
    dropInf <- function(x, h) {
        if (any(isInf <- h >= 1)) {
            warning(gettextf("not plotting observations with leverage one:\n  %s", 
                paste(which(isInf), collapse = ", ")), call. = FALSE, 
                domain = NA)
            x[isInf] <- NaN
        }
        x
    }
    show <- rep(FALSE, 6)
    show[which] <- TRUE
    class(object$residuals) <- NULL
    if (is.null(data)) {
        plot.data <- ggplot2::fortify(object)
    }
    else {
        plot.data <- ggplot2::fortify(object, data = data)
    }
    n <- nrow(plot.data)
    plot.data$.index <- 1:n
    plot.data$.label <- rownames(plot.data)
    is_glm <- inherits(object, "glm")
    r <- residuals(object)
    w <- weights(object)
    if (any(show[2L:6L])) {
        s <- if (inherits(object, "rlm")) {
            object$s
        }
        else if (is_glm) {
            sqrt(summary(object)$dispersion)
        }
        else {
            sqrt(stats::deviance(object)/stats::df.residual(object))
        }
        hii <- stats::lm.influence(object, do.coef = FALSE)$hat
        r.hat <- range(hii, na.rm = TRUE)
        is_const_lev <- all(r.hat == 0) || all(diff(r.hat) < 
            1e-10 * mean(hii, na.rm = TRUE))
        fs <- dplyr::select_if(plot.data, function(x) is.character(x) | 
            is.factor(x))
        fs[[".label"]] <- NULL
        if (is_const_lev & ncol(fs) > 0) {
            plot.data$.nf <- stringr::str_wrap(interaction(fs, 
                sep = ":"), width = 10)
        }
        if (any(show[2L:3L])) {
            plot.data$.wresid <- if (is.null(w)) {
                r
            }
            else {
                sqrt(w) * r
            }
            plot.data$.wstdresid <- plot.data$.wresid/(s * sqrt(1 - 
                hii))
        }
        if (show[2L]) {
            ylim <- range(plot.data$.wstdresid, na.rm = TRUE)
            ylim[2L] <- ylim[2L] + diff(ylim) * 0.075
            qn <- stats::qqnorm(plot.data$.wstdresid, ylim = ylim, 
                plot.it = FALSE)
            plot.data$.qqx <- qn$x
            plot.data$.qqy <- qn$y
        }
    }
    label.fitted <- ifelse(is_glm, "Predicted values", "Fitted values")
    label.y23 <- ifelse(is_glm, "Std. deviance resid.", "Standardized residuals")
    if (is.logical(shape) && !shape) {
        if (missing(label)) {
            label <- TRUE
        }
        if (missing(label.n)) {
            label.n <- nrow(plot.data)
        }
    }
    plot.data <- flatten(plot.data)
    if (label.n > 0L) {
        if (show[1L]) {
            r.data <- dplyr::arrange(plot.data, dplyr::desc(abs(.resid)))
            r.data <- utils::head(r.data, label.n)
        }
        if (".wresid" %in% colnames(plot.data)) {
            wr.data <- dplyr::arrange(plot.data, dplyr::desc(abs(.wresid)))
            wr.data <- utils::head(wr.data, label.n)
        }
        if (any(show[4L:6L])) {
            cd.data <- dplyr::arrange(plot.data, dplyr::desc(abs(.cooksd)))
            cd.data <- utils::head(cd.data, label.n)
        }
    }
    .smooth <- function(x, y) {
        stats::lowess(x, y, f = 2/3, iter = 3)
    }
    .decorate.label <- function(p, data) {
        if (label & label.n > 0) {
            p <- plot_label(p = p, data = data, label = label, 
                label.label = label.label, label.colour = label.colour, 
                label.alpha = label.alpha, label.size = label.size, 
                label.angle = label.angle, label.family = label.family, 
                label.fontface = label.fontface, label.lineheight = label.lineheight, 
                label.hjust = label.hjust, label.vjust = label.vjust, 
                label.repel = label.repel)
        }
        p
    }
    .decorate.plot <- function(p, xlab = NULL, ylab = NULL, title = NULL) {
        p + ggplot2::xlab(xlab) + ggplot2::ylab(ylab) + ggplot2::ggtitle(title)
    }
    smoother_m <- ggplot2::aes_string(x = "x", y = "y")
    if (show[1L]) {
        t1 <- "Residuals vs Fitted"
        mapping <- ggplot2::aes_string(x = ".fitted", y = ".resid")
        smoother <- .smooth(plot.data$.fitted, plot.data$.resid)
        smoother <- as.data.frame(smoother)
        p1 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p1 <- p1 + geom_factory(geom_point, plot.data, colour = colour, 
                size = size, linetype = linetype, alpha = alpha, 
                fill = fill, shape = shape)
        }
        p1 <- p1 + ggplot2::geom_line(data = smoother, mapping = smoother_m, 
            colour = smooth.colour, linetype = smooth.linetype) + 
            ggplot2::geom_hline(yintercept = 0L, linetype = ad.linetype, 
                size = ad.size, colour = ad.colour)
        p1 <- .decorate.label(p1, r.data)
        p1 <- .decorate.plot(p1, xlab = label.fitted, ylab = "Residuals", 
            title = t1)
    }
    if (show[2L]) {
        t2 <- "Normal Q-Q"
        qprobs <- c(0.25, 0.75)
        qy <- stats::quantile(plot.data$.wstdresid, probs = qprobs, 
            names = FALSE, type = 7, na.rm = TRUE)
        qx <- stats::qnorm(qprobs)
        slope <- diff(qy)/diff(qx)
        int <- qy[1L] - slope * qx[1L]
        mapping <- ggplot2::aes_string(x = ".qqx", y = ".qqy")
        p2 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p2 <- p2 + geom_factory(geom_point, plot.data, colour = colour, 
                size = size, linetype = linetype, alpha = alpha, 
                fill = fill, shape = shape)
        }
        p2 <- p2 + ggplot2::geom_abline(intercept = int, slope = slope, 
            linetype = ad.linetype, size = ad.size, colour = ad.colour)
        p2 <- .decorate.label(p2, wr.data)
        p2 <- .decorate.plot(p2, xlab = "Theoretical Quantiles", 
            ylab = label.y23, title = t2)
    }
    if (show[3L]) {
        t3 <- "Scale-Location"
        mapping <- ggplot2::aes_string(x = ".fitted", y = "sqrt(abs(.wstdresid))")
        smoother <- .smooth(plot.data$.fitted, sqrt(abs(plot.data$.wstdresid)))
        smoother <- as.data.frame(smoother)
        p3 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p3 <- p3 + geom_factory(geom_point, plot.data, colour = colour, 
                size = size, linetype = linetype, alpha = alpha, 
                fill = fill, shape = shape)
        }
        p3 <- p3 + ggplot2::geom_line(data = smoother, mapping = smoother_m, 
            colour = smooth.colour, linetype = smooth.linetype)
        p3 <- .decorate.label(p3, wr.data)
        label.y3 <- ifelse(is_glm, expression(sqrt(abs(`Std. deviance resid.`))), 
            expression(sqrt(abs(`Standardized residuals`))))
        p3 <- .decorate.plot(p3, xlab = label.fitted, ylab = label.y3, 
            title = t3)
    }
    if (show[4L]) {
        t4 <- "Cook's distance"
        mapping <- ggplot2::aes_string(x = ".index", y = ".cooksd", 
            ymin = 0, ymax = ".cooksd")
        p4 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p4 <- p4 + geom_factory(geom_linerange, plot.data, 
                colour = colour, size = size, linetype = linetype, 
                alpha = alpha, fill = fill, shape = shape)
        }
        p4 <- .decorate.label(p4, cd.data)
        p4 <- .decorate.plot(p4, xlab = "Obs. Number", ylab = "Cook's distance", 
            title = t4)
    }
    if (show[5L]) {
        if (is_const_lev & ncol(fs) > 0) {
            t5 <- "Constant Leverage:\nResiduals vs Factor Levels"
            mapping <- ggplot2::aes_string(x = ".nf", y = ".stdresid")
            p5 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
            if (!is.logical(shape) || shape) {
                p5 <- p5 + geom_factory(geom_point, plot.data, 
                  colour = colour, size = size, linetype = linetype, 
                  alpha = alpha, fill = fill, shape = shape)
            }
            p5 <- p5 + ggplot2::geom_hline(yintercept = 0L, linetype = ad.linetype, 
                size = ad.size, colour = ad.colour) + ggplot2::expand_limits(x = 0)
            p5 <- .decorate.label(p5, cd.data)
            label.y5 <- ifelse(is_glm, "Std. Pearson resid.", 
                "Standardized Residuals")
            p5 <- .decorate.plot(p5, xlab = "Factor Level Combination", 
                ylab = label.y5, title = t5)
        }
        else {
            t5 <- "Residuals vs Leverage"
            mapping <- ggplot2::aes_string(x = ".hat", y = ".stdresid")
            smoother <- .smooth(plot.data$.hat, plot.data$.stdresid)
            smoother <- as.data.frame(smoother)
            p5 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
            if (!is.logical(shape) || shape) {
                p5 <- p5 + geom_factory(geom_point, plot.data, 
                  colour = colour, size = size, linetype = linetype, 
                  alpha = alpha, fill = fill, shape = shape)
            }
            p5 <- p5 + ggplot2::geom_line(data = smoother, mapping = smoother_m, 
                colour = smooth.colour, linetype = smooth.linetype) + 
                ggplot2::geom_hline(yintercept = 0L, linetype = ad.linetype, 
                  size = ad.size, colour = ad.colour) + ggplot2::expand_limits(x = 0)
            p5 <- .decorate.label(p5, cd.data)
            label.y5 <- ifelse(is_glm, "Std. Pearson resid.", 
                "Standardized Residuals")
            p5 <- .decorate.plot(p5, xlab = "Leverage", ylab = label.y5, 
                title = t5)
        }
    }
    if (show[6L]) {
        t6 <- "Cook's dist vs Leverage"
        mapping <- ggplot2::aes_string(x = ".hat", y = ".cooksd")
        smoother <- .smooth(plot.data$.hat, plot.data$.cooksd)
        smoother <- as.data.frame(smoother)
        p6 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p6 <- p6 + geom_factory(geom_point, plot.data, colour = colour, 
                size = size, linetype = linetype, alpha = alpha, 
                fill = fill, shape = shape)
        }
        p6 <- p6 + ggplot2::geom_line(data = smoother, mapping = smoother_m, 
            colour = smooth.colour, linetype = smooth.linetype) + 
            ggplot2::expand_limits(x = 0, y = 0)
        p6 <- .decorate.label(p6, cd.data)
        p6 <- .decorate.plot(p6, xlab = "Leverage", ylab = "Cook's distance", 
            title = t6)
        g <- dropInf(hii/(1 - hii), hii)
        p <- length(stats::coef(object))
        bval <- pretty(sqrt(p * plot.data$.cooksd/g), 5)
        for (i in seq_along(bval)) {
            bi2 <- bval[i]^2
            p6 <- p6 + ggplot2::geom_abline(intercept = 0, slope = bi2, 
                linetype = ad.linetype, size = ad.size, colour = ad.colour)
        }
    }
    if (is.null(ncol)) {
        ncol <- 0
    }
    if (is.null(nrow)) {
        nrow <- 0
    }
    plot.list <- list(p1, p2, p3, p4, p5, p6)[which]
    new("ggmultiplot", plots = plot.list, nrow = nrow, ncol = ncol)
}
<bytecode: 0x609dfc00dcb8>
<environment: namespace:ggfortify>

sloop::s3_get_method(augment.lm)

function (x, data = model.frame(x), newdata = NULL, se_fit = FALSE, 
    interval = c("none", "confidence", "prediction"), conf.level = 0.95, 
    ...) 
{
    warn_on_subclass(x, "augment")
    check_ellipses("level", "augment", "lm", ...)
    interval <- match.arg(interval)
    df <- augment_newdata(x, data, newdata, se_fit, interval, 
        level = conf.level)
    if (is.null(newdata)) {
        tryCatch({
            infl <- influence(x, do.coef = FALSE)
            df <- add_hat_sigma_cols(df, x, infl)
        }, error = data_error)
    }
    df
}
<bytecode: 0x609dfba87d00>
<environment: namespace:broom>

Question

Solution

sloop::s3_methods_class("prcomp") |>
  gt::gt() |>
  gt::tab_header("Generics of class prcomp")

Generics of class prcomp
generic	class	visible	source
augment	prcomp	FALSE	registered S3method
autoplot	prcomp	FALSE	registered S3method
biplot	prcomp	FALSE	registered S3method
fortify	prcomp	FALSE	registered S3method
plot	prcomp	FALSE	registered S3method
predict	prcomp	FALSE	registered S3method
print	prcomp	FALSE	registered S3method
summary	prcomp	FALSE	registered S3method
tidy	prcomp	FALSE	registered S3method

Question

Solution

sloop::s3_dispatch(autoplot(lm0))

=> autoplot.lm
 * autoplot.default

Question

Solution

data(UCBAdmissions)
class(UCBAdmissions)

[1] "table"

otype(UCBAdmissions)

[1] "S3"

Question

Solution

sloop::s3_class(UCBAdmissions)

[1] "table"

# sloop::s3_get_method(autoplot.table)
# sloop::s3_get_method(augment.table)

• OO in Advanced R Programming 1st Edition
• S3 classes Advanced R 2nd Edition

Programming with dplyr and ggplot2

We first aim at programming a function that takes as input a dataframe df, a column name col, and that, depending on the type of the column denoted by col, plots a histogram (for numerical column), a barplot (for factors), or raise an error of the column is neither categorical, nor numerical.

The function should return a ggplot object.

Here is a first attempt.

Let us first build a toy tibble.

tb <- tibble( 
  col_num = rnorm(100), 
  col_fac = as_factor(sample(letters, 100, replace = T)), 
  col_ts = Sys.time() + duration(sample(1:20, 100, replace=T),units="days")
) 

tb |> 
  head()

# A tibble: 6 × 3
  col_num col_fac col_ts             
    <dbl> <fct>   <dttm>             
1  -0.425 z       2025-04-04 00:11:08
2   0.129 l       2025-03-23 23:11:08
3  -0.593 i       2025-03-23 23:11:08
4   0.554 s       2025-04-02 00:11:08
5  -1.21  e       2025-03-29 23:11:08
6   0.273 k       2025-03-25 23:11:08

Solution

gg_obj <-  function(df, col){
  
  vct <- df[[col]]
  tp <- class(vct)

  if (tp != "numeric" & tp !="factor") {
    stop(paste0(col, " is of wrong type!"))
  }

  p <- ggplot(df) + 
    aes(x=.data[[col]])

  if (tp=="numeric") {
    p <- p + geom_histogram()
  } else {
    p <- p + geom_bar()
  }

  p  
}

1

List component is accessed by name. col is a string.

2

.data is a pronoun for the dataframe component of the ggplot object

(
  gg_obj(tb, "col_num") +
    labs(
      title= "Histogram", 
      subtitle= "Numerical column")
) + (
  gg_obj(tb, "col_fac") +
    labs(
      title= "Barplot",
      subtitle= "Factor column"
    )
)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Question

• Pass more optional arguments to geom_... (use ellipsis ...)
• Avoid quoting the column name

gg_obj_2 <-  function(df, col, ...){
  params <- list(...)
  vct <- pull(df, {{col}})
  tp <- class(vct)[1]

  if (tp != "numeric" & tp !="factor") {
    stop("column is of wrong type!")
    return
  }

  p <- ggplot(df) + 
    aes(x={{col}})

  if (tp=="numeric") {
    p <- p + geom_histogram(...)
  } else {
    p <- p + geom_bar(...)
  }

  p  
}

1

standard technique for programming with tidy-evaluation: embracing.

2

using the ellipsis

(
  gg_obj(tb, "col_num") +
    labs(
      title= "Histogram", 
      subtitle= "Numerical column")
) +
(
  gg_obj_2(tb, col_num, fill="red") +
    labs(
      title= "Histogram, tuning color", 
      subtitle= "Numerical column")
)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Question

Question

How could you add a geom_point() layer to each element of the following list?

plots <- list(
  ggplot(mpg, aes(displ, hwy)),
  ggplot(diamonds, aes(carat, price)),
  ggplot(faithfuld, aes(waiting, eruptions, size = density))
)

From R Advanced Programming

Solution

plots |> 
  map(\(x) x + geom_point()) |>
  patchwork::wrap_plots()

Inside lm()

Question

In classes like lm, prcomp, … we have a member called call. What does it represent? How is it constructed?

First, read the code of lm.

> lm 
function (formula, data, subset, weights, na.action, method = "qr", 
    model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, 
    contrasts = NULL, offset, ...) 
{
    ret.x <- x
    ret.y <- y
    cl <- match.call()
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("formula", "data", "subset", "weights", "na.action", 
        "offset"), names(mf), 0L)
    mf <- mf[c(1L, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1L]] <- quote(stats::model.frame)
    mf <- eval(mf, parent.frame())
    if (method == "model.frame") 
        return(mf)
    else if (method != "qr") 
        warning(gettextf("method = '%s' is not supported. Using 'qr'", 
            method), domain = NA)
    mt <- attr(mf, "terms")
    y <- model.response(mf, "numeric")
    w <- as.vector(model.weights(mf))
    if (!is.null(w) && !is.numeric(w)) 
        stop("'weights' must be a numeric vector")
    offset <- model.offset(mf)
    mlm <- is.matrix(y)
    ny <- if (mlm) 
        nrow(y)
    else length(y)
    if (!is.null(offset)) {
        if (!mlm) 
            offset <- as.vector(offset)
        if (NROW(offset) != ny) 
            stop(gettextf("number of offsets is %d, should equal %d (number of observations)", 
                NROW(offset), ny), domain = NA)
    }
    if (is.empty.model(mt)) {
        x <- NULL
        z <- list(coefficients = if (mlm) matrix(NA_real_, 0, 
            ncol(y)) else numeric(), 
                  residuals = y, 
                  fitted.values = 0 * y, 
                  weights = w, 
                  rank = 0L, 
                  df.residual = if (!is.null(w)) sum(w != 0) else ny
              )
        if (!is.null(offset)) {
            z$fitted.values <- offset
            z$residuals <- y - offset
        }
    }
    else {
        x <- model.matrix(mt, mf, contrasts)
        z <- if (is.null(w)) 
            lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)
        else lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok, 
            ...)
    }
    class(z) <- c(if (mlm) "mlm", "lm")
    z$na.action <- attr(mf, "na.action")
    z$offset <- offset
    z$contrasts <- attr(x, "contrasts")
    z$xlevels <- .getXlevels(mt, mf)
    z$call <- cl
    z$terms <- mt
    if (model) 
        z$model <- mf
    if (ret.x) 
        z$x <- x
    if (ret.y) 
        z$y <- y
    if (!qr) 
        z$qr <- NULL
    z
}
<bytecode: 0x55564224e930>
<environment: namespace:stats>

Question

Have a look at function match.call()

e <- match.call(get, call("get", "abc", i = FALSE, p = 3))

is_expression(e)

[1] TRUE

is_call(e)

[1] TRUE

lobstr::ast(match.call(get, call("get", "abc", i = FALSE, p = 3)))

█─match.call 
├─get 
└─█─call 
  ├─"get" 
  ├─"abc" 
  ├─i = FALSE 
  └─p = 3 

Let us explore the next toy example.

fun <- function(x, lower = 0, upper = 1) {
  structure((x - lower) / (upper - lower), CALL = match.call())
}


w <- fun(4 * atan(1), u = pi)

type_of(w)

Warning: `type_of()` is deprecated as of rlang 0.4.0.
Please use `typeof()` or your own version instead.
This warning is displayed once every 8 hours.

[1] "double"

w

[1] 1
attr(,"CALL")
fun(x = 4 * atan(1), upper = pi)

str(w)

 num 1
 - attr(*, "CALL")= language fun(x = 4 * atan(1), upper = pi)

x <- attr(w, 'CALL')

typeof(x)

[1] "language"

is_call(x)

[1] TRUE

eval(x)  

[1] 1
attr(,"CALL")
fun(x = 4 * atan(1), upper = pi)

# try also eval_tidy(), eval_bare

Data masking and environments

Question

Solution

Tidy evaluation

Question

What is quasi-quotation?

Keep the rlang cheatsheet around.

Question

Explain the difference between an expression and a quosure

Question

References

Programming with ggplot